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Computing  Locally-Mass-Conservative 
Fluxes  from  Multi-dimensional  Finite 
Element  Flow  Simulations 


by  Jing-Ru  C.  Cheng,  Hwai-Ping  Cheng,  Matthew  W.  Farthing,  and  Christopher  E.  Kees 


PURPOSE:  This  technical  note  details  how  conservative  normal  fluxes  through  element  edges 
and  faces  are  computed  in  two-  (2-D)  and  three-dimensional  (3-D)  spaces  in  a  flux-calculation 
post-processor  developed  in  the  System-Wide  Water  Resources  Program  (SWWRP). 

BACKGROUND:  Conserving  local  mass  in  the  finite  volume  (FV)  sense,  where  the  sum  of 
face  fluxes  equal  to  the  rate  of  change  of  storage  within  each  element/cell,  is  essential  in  comput¬ 
ing  water  flow  and  contaminant  transport.  Although  the  continuous  Galerkin  finite  element  (FE) 
method  does  not  yield  locally  conservative  flux  approximation  directly,  Berger  and  Howington 
showed  that,  by  remaining  consistent  with  the  discrete  approximation  given  by  the  FE  statement, 
the  resulting  flux  estimates  will  preserve  mass  balance  [1].  Passing  conservative  water  flux 
through  each  element  edge/face  from  flow  models  to  transport  models  is  critical  for  accurate 
simulation  and  analysis.  At  the  U.S.  Army  Engineer  Research  and  Development  Center  (ERDC), 
most  water  flow  models  employ  the  FE  method,  while  many  contaminant  transport  models  use 
the  FV  method.  The  computation  of  locally-conservative  water  flux  through  each  elemental 
edge/face  has  thus  become  necessary  for  passing  FE-based  water  flux  to  FV-based  contaminant 
transport  models.  This  technical  note  describes  how,  from  multi-dimensional  FE-based  flow 
simulations,  we  computed  locally-mass-conservative  fluxes  to  hand-off  to  FV-based  models  and 
to  eliminate  apparent  flux  jumps  on  element  boundaries  for  particle  tracking.  The  computed  con¬ 
servative  flux  can  also  be  used  to  set  boundary  conditions  for  inset  models  when  desired. 

METHODOLOGY:  We  employed  the  two  flux  post-processing  techniques  focused  in  [2].  For 
convenience,  we  refer  to  the  first  technique  as  the  local  approach,  and  the  second  one  the  global 
approach  in  this  technical  note.  Both  were  originally  developed  by  Larson  and  Niklasson  [3]. 
Larson  and  Niklasson’s  global  approach  is  mathematically  equivalent  to  Sun  and  Wheeler’s 
global  approach  [4].  For  both  the  local  and  the  global  approaches,  the  required  input  data  include 

(1)  the  element  residuals  associated  with  the  respective  computed  solution,  meaning  they  are 
computed  by  constructing  matrix  equations  at  the  element  level  with  flux-type  boundary 
conditions  incorporated  (note  that  each  element  has  residuals  matching  its  element  nodes 
to  satisfy  local  conservation); 

(2)  boundary  face  information,  including 

a.  an  indicator  to  show  whether  this  boundary  face  is  specified  with  a  flux-type  boun¬ 
dary  condition: 

•  a  no-flow  boundary  face  is  specified  with  a  zero-flux  boundary  condition, 
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•  a  flow-through  (open)  boundary  face  specified  with  a  flux-type  boundary  condi¬ 
tion  will  have  the  normal  flux  match  the  specified  value, 

•  a  flow-through  (open)  boundary  face  without  a  flux-type  boundary  condition 
applied  will  have  the  normal  flux  computed. 

b.  the  specified  normal  flux  for  this  boundary  face  (set  to  zero  if  the  normal  flux  through 
this  boundary  face  is  not  to  be  computed) 

(3)  the  estimated  fluxes  associated  with  element  faces; 

(4)  the  element  connectivity  infonnation; 

(5)  nodal  coordinates. 

From  the  input  data  (3)  -  (5)  above,  we  can  derive  the  following  information  for  conservative 
nonnal  flux  computation: 

(1)  node-element  connectivity, 

(2)  edge-element  (2-D)  or  face-element  (3-D)  connectivity, 

(3)  edge-node  (2-D)  or  face-node  (3-D)  connectivity, 

(4)  edge  length  (2-D)  or  face  area  (3-D). 

It  is  noted  that  a  no-flow  boundary  face,  by  definition,  is  associated  with  a  zero-flux  boundary 
condition,  and  thus  the  normal  flux  through  this  boundary  face  is  zero.  On  the  other  hand,  a 
boundary  face  is  flow-through  if  it  is  not  associated  with  a  no-flow  boundary  condition.  A  flow¬ 
through  boundary  face  can  be  assigned  a  non-zero  flux-type  boundary  condition,  or  its  associated 
boundary  flux  is  computed  based  on  the  solution  of  the  continuity  equation.  This  concept  is 
applied  when  we  compute  conservative  fluxes  through  element  faces  associated  with  boundary 
nodes. 

Figure  1  depicts  the  flow  chart  for  computing  conservative  normal  fluxes  in  the  post-processor: 
reading  input  data  in  Step  1,  followed  by  generating  connectivity  information,  face  length/area, 
and  boundary  flux  indicator  in  Step  2,  then  the  computation  of  an  artificial,  element-based  quan¬ 
tity  used  for  nonnal  flux  correction  in  Step  3,  and  finally  in  Step  4  the  computation  of  conserva¬ 
tive  normal  fluxes  through  all  element  faces.  In  the  following,  we  discuss  the  computation 
involved  in  Steps  3  and  4  with  the  local  and  the  global  approaches. 
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Figure  1 .  Flow  chart  for  computing  conservative  normal  flux.  “Face”  is  used  in  the  chart  to  represent 
2-D  edge  or  3-D  face. 


The  Local  Approach.  In  the  local  approach,  the  conservative  sub-edge  (in  2-D)  or  sub-face 
(in  3-D)1  normal  fluxes  associated  with  element  faces  around  each  node  are  computed  on  a  node- 
by-node  basis.  The  conservative  normal  flux  of  an  element  edge/face  is  computed  then  by  com¬ 
bining  all  of  the  associated  sub-edge  or  sub-face  normal  fluxes. 

Scenario  1:  For  element  edges  associated  with  an  internal  node.  As  shown  in  Figure  2  below, 
there  are  6  elements  connected  at  Node  1,  and  there  are  six  sub-edge  flows  satisfying  the  discre¬ 
tized  governing  equation  of  continuity,  e.g.,  Richards’  equation.  In  Figure  2(b),  the  three  resi¬ 
duals,  each  associated  with  a  node,  of  each  element  are  marked.  These  residuals  represent  the  net 
flow  entering  or  leaving  the  element,  and  they  can  be  computed  by  substituting  the  numerical 


1  A  sub-edge  or  a  sub-face  of  a  node  is  an  elemental  edge/face  containing  the  node.  A  sub-edge  has  a  length  equal  to 
half  of  the  element  edge  length,  while  a  sub-face  has  an  area  equal  to  the  element  face  area  divided  by  the  number  of 
face  nodes.  It  is  called  a  sub-edge  or  a  sub-face  because  the  conservative  normal  flux  computed  for  a  sub-edge  or 
sub-face  is  associated  with  only  a  portion  of  edge  length  or  face  area.  For  example,  in  Figure  2,  there  are  six  element 
edges  around  Node  1: 1 is,  I  is,  I  is,  I  is,  I  is,  1 1-7 ■  While  each  of  these  six  edges  is  a  sub-edge  of  Node  1,  each  sub¬ 
edge  has  a  length  only  half  of  the  corresponding  element  edge  length. 
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Figure  2.  The  six  triangular  elements  associated  with  an  internal  node  (i.e.,  Node  1)  and  the  six 

conservative  sub-edge  flows  around  Node  1 :  (a)  elements  are  connected,  (b)  elements  are 
separated. 


solution  of  the  continuity  equation  back  into  the  element  matrix  equation  constructed  in  the  FE 
integration.  For  example,  if  we  solve  the  Richards’  equation  for  subsurface  flow  using  the  Galer- 
kin  FE  method,  the  element  residual  at  local  node  j  can  be  defined  as 


R‘: 


where 

Rj  = 

a  = 

N‘  = 
F  = 
hj  = 
V  = 

K  = 

q  = 

h  = 
z  = 
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=  -\n-K-(yh  +  Wz)NeidTe 
re 


f  NfFN'dn, 
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J  '  J  e 

Qe 

dt 

J  V  1  )  \  j )  e 

ae 

the  element  residual  at  local  node  j 

the  domain  of  element  e 

the  i-th  local  base  function  of  element  e 

water  capacity  associated  with  element  e 
pressure  head  at  local  node  j 

the  del  operator 

hydraulic  conductivity  tensor  associated  with  element  e 
source/sink  in  element  e 

pressure  head 
gravity  head 

the  boundary  of  element  e 


(1) 
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Here,  Rf  represents  the  net  outgoing  flow  associated  with  local  node  j  when  element  e  is  consi¬ 
dered.  Therefore  in  Figure  2(b),  R{4,  R-i,  and  Rf  are  the  three  element  residuals  of  Element  (A); 
Ri  ,  i?3  ,  and  R4  are  the  three  element  residuals  of  Element  ( B ),  and  so  on.  Because  Node  1  is  an 
internal  node  and  there  are  no  sources  or  sinks,  the  sum  of  all  the  six  element  residuals  associated 
with  Node  1  is  zero,  i.e., 

/?,  =  Rf  +  Rf  +  Rf  +  Rf  +  Rf  +Rf  =  0  (2) 


where 

Ri  =  total  element  residual  associated  with  Node  i 
Rf  =  element  residual  associated  with  Node  i  coming  from  element  (E) 


and,  Rf  can  be  expressed  as 

p  A  rmc  rmc  J/mc  ^  AB  ]/ ttlc  ^ FA 

~  Jl.AB  _  J\,FA  ~  '  1,AB  '  V\ ,FA  '  ~ 


(3) 


where 

=  nonnal  sub-edge  flow  from  Element  (Ef)  to  Element  (E2),  which  conserves  mass 
locally  around  Node  i 

VffE2  =  normal  sub-edge  flux  from  Element  (£/)  to  Element  (E2),  which  conserves  mass 
locally  around  Node  i 

/££  =  edge  length  associated  with  the  interface  between  Elements  (£/)  and  (E2) 


If  Ei  .in  represents  a  given  estimated  normal  sub-edge  flux,  we  can  relate  Vff  to  V\  jb  with  a 
correction  term  as 


va,  =  Vui  +  A  V,M  =  Vt  AB  +(uu  -  UI  A ) 


(4) 


where  U\^  and  U\  n  are  artificial,  element-based  quantities  introduced  to  help  calculate  the  cor¬ 
rection  of  normal  sub-edge  flux. 

Therefore,  Equation  3  can  be  written  as: 


R?  = 


y,.FA+(uu-u,F) 


lFA 

2 


(5) 


Likewise,  Rf ,  Rf ,  Rf ,  Rf ,  Rf  can  be  expressed  as 


Rf  = 


V+Kc-t/,.*)}y 


(6) 
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Ken  +  (K,n  ~ K,e )]~-[Kbc  +  (K,c-Ukb) 


‘bc 

2 


RT  = 


Kcd+{Uw-u,.c) 


L CD 

2 


Rf  = 


i 


KEF+{K,F-uhE)  vUDE+(uhE-uhD) 


[de_ 
J  2 


R  = 


\  FA+(UU-UIF) 


lFA 

2 


y  +(tj  —  TJ 
v\iEF^yu\iF  L'/l,£/J  2 


(7) 

(8) 

(9) 

(10) 


It  must  be  noted  that  every  equation  from  Equation  5  through  Equation  10  can  be  represented  by 
the  other  five  equations,  i.e.,  there  are  only  five  independent  equations.  To  solve  the  six 
unknowns,  we  enforce  UlA  =0  and  use  it  to  replace  one  of  the  six  equations,  say  Equation  5. 

This  will  not  affect  the  resulting  flux  correction  mathematically  because  Equations  6  through 
Equation  10  correspond  to  a  Neumann  problem,  as  mentioned  in  Kees  et  al.  [2].  It  is  because  the 
flux  correction  is  the  jump  of  two  adjacent  U’s,  rather  than  the  U  values.  Setting  U/j  to  any 
value  will  still  produce  the  same  jump  values  as  long  as  the  element  residual  equations  are  solved 
accurately. 

After  solving  Equation  6  through  Equation  10,  the  conservative  normal  sub-edge  fluxes  can  be 
computed  with 


(11a) 

Kbc=V,.bc+{Uu-Uu ,) 

(lib) 

ymc  —V  -\-{u  —  U  ^ 

V  \,CD  V  \,CD  ^  V  1  >D  ^  1,C  ) 

(He) 

vr»E=vim+(u,j,-utJ>) 

(Hd) 

v;%=rtl:F+(u^F-ulE) 

(lie) 

Kfa=Kfa+{0-Uff) 

(Ilf) 

Scenario  2:  For  element  edges  associated  with  a  boundary  node.  A  boundary  edge  can  be  of 
either  flow-through  or  no-flow.  Each  flow-through  boundary  edge  may  or  may  not  be  assigned  a 
specified  flux-type  boundary  condition  in  numerical  simulation.  When  it  is  assigned  with  a  speci¬ 
fied  flux-type  boundary  condition,  its  normal  flux  must  be  identical  to  the  specified  flux.  If  it  is 
not  specified  with  a  flux-type  boundary  condition,  its  nonnal  flux  can  be  computed  based  on  the 
solution  of  the  continuity  equation.  A  no-flow  boundary  edge,  on  the  other  hand,  is  associated 
with  a  zero  nonnal  flux  by  definition.  In  other  words,  it  is  specified  with  a  zero-flux  boundary 
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condition  even  though  in  most  numerical  models  a  boundary  edge  is  usually  defaulted  to  no-flow 
type  when  it  is  not  specified  with  any  type  of  boundary  condition. 

As  shown  in  Figure  3,  there  are  three  elements  around  Node  1,  which  is  a  boundary  node,  and 
Edges  1-2  and  1-5  are  the  two  boundary  edges  adjacent  to  Node  1.  In  this  case,  there  are  four 
normal  sub-edge  fluxes  to  be  computed  to  satisfy  the  discretized  continuity  equation. 


Figure  3.  The  three  triangular  elements  associated  with  a  boundary  node  (i.e.,  Node  1)  and  the  four 
conservative  sub-edge  flows  around  Node  1 :  (a)  elements  are  connected,  (b)  elements  are 
separated;  Among  the  four  sub-edge  flows,  two  are  associated  with  boundary  edges. 


Without  the  flux -type  boundary  conditions  applied,  the  three  element  residual  equations  asso¬ 
ciated  with  Node  1  are; 


n  A  f'ne  _  fine  _  ymc  ^  AB  ymc  ^ ' ZA 

Al  ~  J\,AB  J\,2A~V\,Ab'  2  V\,ZA  ‘  2 


Xa+l^u-^u)]- y 


TfB  _  fine  _  fine  _  ymc  ^  BC  _  y  me  l AB 

Al  ~  J\,BC  J\,AB~Vl,BC'  2  V1,AB  ' 


l 


Kbc + (<v  - )  I  ~  -  k,. + (w, -ut,A ) 


1  ^  AB 

J  2 


nC  _  rmc  fmc  -trine  ^CX Trine  ^ BC 

“  J\,cx  J \,bc  yi,cx  2  '  2 


,cx 


J  2 


2 


(12) 


(13) 


(14) 
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where  fmZA  and  fmfx  are  the  sub-edge  flow  through  Boundary  Edges  1-2  and  1-5,  respectively; 
lZA  and  lcx  are  the  edge  lengths  associated  with  Boundary  Edges  1-2  and  1-5,  respectively;  Ux  z 
and  Ux  x  are  the  U  values  associated  imaginary  elements  outside  of  the  computational  domain. 


Without  knowing  the  relationship  between  Ux  x  and  Ulz,  we  can  simply  assume  and  enforce 
UIX  =  Ulz  =C  in  Equation  12  and  Equation  14  to  help  close  the  computational  system,  where  C 

is  a  constant.  The  correction  of  normal  sub-edge  flux  (i.e.,  the  jump  of  El  between  connected 
elements)  is  uniquely  determined  as  long  as  Ux  x  and  Ui  z  are  set  to  the  same  value,  which  is 

reasonable  because  they  represent  outside  quantities  around  Node  1.  By  setting  C  =  0, 
Equations  12  through  14  thus  become 


Rf  =  fZ*  -  fZCA  =  [K,ab + {ulB  [vlZA  +(ulA-o) 


lZA 

2 


»B  rmc  _  rmc  fy  ,  ( TJ  _TJ  )j  .  [§£.  -  [v  +  ( TJ  —TJ  V  ^AB 
1V\  J\,BC  J\,AB  \_V  \,BC  ^  V  1>C  J J  9  \,AB  ^  \,B  U1  ,A ) _ 


K=f,7x-fw=[V,.cx+(°-U>.c) 


lcx 

2 


J  2 


k,bc+(  u^-u^y-f 


(15) 

(16) 

(17) 


In  solving  Equation  15  through  Equation  17,  we  consider  three  possible  situations  concerning  the 
specified  flux-type  boundary  condition:  (1)  when  all  boundary  edges  are  assigned;  (2)  when 
some,  but  not  all,  boundary  edges  are  assigned;  (3)  when  no  boundary  edge  is  assigned. 

Situation  1.  All  boundary  edges  are  assigned  flux-type  boundary  conditions:  When  flux-type 
boundary  conditions  have  been  taken  into  account  in  computing  the  element  residuals,  f"ZA  and 

fx'cx  will  no  longer  appear  in  the  element  residual  equations.  If  qzA  and  qcx  represent  the  out¬ 
ward  boundary  flux  through  Edges  1-2  and  1-5,  respectively,  Equations  15  through  17  become: 


pA'bf  =  pA  _^zajJza_ 

1  r\ 


lab 

2 


R°’bf  =  R"  =  [Kbc  +  (uUc  [Vuab  + (Uu,  -  UhA  ) 


i AB 

/J  2 


pC,bf  _  rC  _  dcx  ‘lex 


r,jc+(u,.c-u,,Bj 


1 BC 
2 


(18) 

(19) 

(20) 


where  RfhJ ,  R'f  h< ,  and  R[  J)J  represent  element  residuals  with  flux-type  boundary  conditions 
taken  into  account. 
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Equation  18  is  linearly  dependent  on  Equation  19  and  Equation  20.  This  is  again  a  Neumann 
problem,  and  we  can  enforce  UXA=0  and  solve  Equation  19  and  Equation  20  for  Ux  B  and  Ut  c . 

The  normal  sub-edge  fluxes  can  be  computed  then  with 


vCza  =  -fe 

v":,=Ka,+{u,.,) 

ymc  —  y  +(u  —  U  1 

v  \,BC  v\ ,5Ct\u1,C  \,B  ) 

t ymc  _ 

'l,CX  —  dcx 


(21a) 
(21b) 
(21c) 
(2  Id) 


Situation  2.  Some,  but  not  all,  boundary  edges  are  assigned  flux-type  boundary  conditions:  In  the 
two-dimensional  case  here,  we  consider  when  only  one  of  the  two  connected  boundary  edges  is 
assigned  a  flux-type  boundary  condition.  Suppose  the  outward  boundary  flux  through  Edge  1-2 
is  qzA,  Equation  15  through  Equation  17  become: 


ram  _  ra  _ 


dzA  '  l ZA 


/ 


AB 

2 


=^=k,c+(^,c  -UhB) 


lbc 

2 


V,AB+(U,t-U,Aj\- ‘-f 


RP  =  R,C  =  [>,„  +  (0  -  C/,  c  )]  •  -y  f,BC  +(U,.C  ~Ult) 


'l_l  ^BC 
/J  2 


(22) 

(23) 

(24) 


We  thus  solve  Equation  22  through  Equation  24  for  U IA,  UXB,  and  Ux  c .  Then  the  normal  sub¬ 
edge  fluxes  can  be  computed  with 


KzA=-qZA  (25a) 

Kt,=vuB+(r,.,-uu)  (25b) 

yttc  =  K.c+(U^-UP  (25c) 

K~x=K,cx+(° -Ut,c)  (25d) 


Situation  3.  None  of  the  boundary  edges  is  assigned  a  flux-type  boundary  condition:  In  this  case, 
we  solve  Equation  15  through  Equation  17  for  UX  A,  UX  B,  and  Ul  c,  and  compute  the  normal 

sub-edge  fluxes  by  using 

r,Z  =  Kz,+(u,,A-o)  (26a) 
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K^  =  vu.+(u,.«~uf 

(26b) 

Vw  =  KJc+{Uia-Uib) 

(26c) 

Kcx=V,.cx+{0-Uia) 

(26d) 

Summary  of  sub-edge  flux  computation.  Suppose  there  are  M  elements  connected  at  a 
node,  the  following  can  be  drawn  from  the  discussion  above  in  the  process  of  constructing  the 
element  residual  equations  around  the  node: 

(1)  Use  element  residuals  from  the  FE  integration  to  construct  the  M  element  residual  equa¬ 
tions,  one  for  each  element. 

(2)  If  the  node  is  an  internal  node,  set  one  U  value  to  zero  and  solve  the  others. 

(3)  If  the  node  is  a  boundary  node,  set  the  “outside”  IT’s  to  zero  (e.g.,  Uxx  =  Ulz  =  0 ). 

(4)  If  the  node  is  a  boundary  node  and  all  the  connected  boundary  edges  are  assigned  flux- 
type  boundary  conditions,  set  one  U  value  to  zero  and  solve  the  others. 

The  constructed  element  residual  equations  are  then  solved  with  a  full-pivoting  direct  solver  to 
minimize  possible  numerical  error.  The  conservative  sub-edge  fluxes  are  computed  using  the 
computed  If  s  as  well  as  the  specified  boundary  fluxes. 

Compute  edge  fluxes  based  on  sub-edge  fluxes.  The  conservative  normal  flux  through 
an  element  edge  can  then  be  computed  by  taking  arithmetic  average  of  its  two  corresponding 
sub-edge  fluxes.  In  Figure  2,  for  example,  the  two  conservative  normal  sub-edge  fluxes  of 
Edge  1-2  are  V™‘AB  and  V"‘A/i  .  The  conservative  normal  flux  of  the  edge  (i.e.,  VAB  )  is  computed  as 

tz  me  .  tz  me 
t ytnc  1  ,AB  2,AB 

V AB  ~  ^ 

Extension  to  three-dimensional  space.  The  extension  of  the  local  approach  from  2-  to 
3-dimensional  space  is  straightforward.  For  easy  visualization,  we  use  hexahedral  elements  for 
demonstration.  By  using  the  summary  of  sub-face  flux  computation  above,  we  construct  the  3-D 
element  residual  equations  for  the  corresponding  scenario  as  follows. 

Scenario  1.  Around  an  internal  node:  As  shown  in  Figure  4,  Node  14  is  an  internal  node  con¬ 
nected  to  eight  hexahedral  elements  and  associated  with  12  sub-face  fluxes.  The  eight  element 
residual  equations  associated  with  Node  14  are 

UH,A=  0  (28) 
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(a)  (b) 


Figure  4.  The  eight  hexahedral  elements  associated  with  an  internal  node  (i.e. ,  Node  14)  and  the 

12  conservative  sub-face  flows  around  Node  14:  (a)  elements  are  connected,  (b)  elements  are 
separated. 


r)  B  _  rmc  rmc  r 

^14  —  J\4,BA  ~  Jl4,FB  J 1 


14  ,BD 


^I4,BA  +(^14,A  ^14, b) 

P 14,BD  +(^14 ,D  ~^14,b) 


A 


BA 


v  +(u  —jj  )~1 .  ^ fb 

y\4,FB  Ttu14,S  ^  \4,F  ) \  ^ 


A 


BD 


(29) 


-)C  _  rmc  rmc  rmc 

M4  —  J \4,DC  ~  J\4,GC  ~  J\4,AC 


A 


DC 


^14,DC  +(^14,C 


V  +(u  -U  )~|  •  ^AC 

yl4,AC^\ul4,C  ^l4,A)j  ^ 


v  +(u  -u  y\ ■  ^ gc 

M4.GC  T\R14,C  ^  14, G  )  J  j 


(30) 


p  D  rmc  rmc  rmc 

\4  ~  J\4,DC  ~  J\4,HD  ~  J\4,BD 


^14 ,DC  +  {P\4,C  ~G\4,D  ) 


A 


DC 


V  +(u  -U  )~|  •  ^HD 

r  \4,HD  T  )u14,i  '-'14,//)J  ^ 


^14, BD  +  {U\4,D  B  ) 


A 


BD 


(31) 


nE  _  rmc  .  rmc  .  r 

K\4  ~  —J\4,FE  +  J\4,EA  ^  J 1 


14,  EG 


E\4'FF  (j~^\4,E 


A 


FE 


+ 


y  +(tj  —  TJ  )~| ,  ^EA 

y\4,EA  T  '-'l4,£)J  ^ 


y  +(jj  _[/  )~|  ,^g_ 

y  14, EG  T \^14,G  '-'l4,£)J  j 


(32) 
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p F  rmc  .  rmc  .  rmc 

K\4  —  J\4,FE  +  J\4,FB  +  J\4,FH 

=  V\4,FE  +  (^14 ,E  ~  ^14, F  ) 

A FE  | 

4  L 

v  +(TJ  U  )1 .  ^FB 

y  \4,FB  ^  \4,B  u14/jJ  j 

(33) 

+  ^14 ,FH  +  {U\4,H  ~  ^14 ,F 

\1  Afh 
>\  4 

p G  rmc  ,  rmc  rmc 

A14  —  —J\4,HG  +  J 14, GC  ~  J\4,EG 

=  ~  V 14, HG  +  (j~^\4,G  ~  ^14 ,H  ) 

a 

1 - 1 

+  [v  +(tj  tj  )] .  4 Gc- 

^  |_  v  14, GC  ^  V  14’c  14, G  )  J  /| 

(34) 

~  ^14 ,EG  +  {U\4 ,G  _  UU'E  ) 

''fr 

1 - 1 

pH  _  rmc  .  rmc  rmc 

A14  —  J\4,HG  +  J\4,HD  ~  J\4,FH 

=  V14  hg  +  (^14 ,G  —  ^14 ,H  ) 

^ HG  + 

4 

F  +(u  u  'll  • 

'  \4,HD  '  14, D  ^  14, H  J  J  ^ 

(35) 

_  V\4,FH  +  14, H  F 

\1  afh 

IJ  4 

where  Aba  denotes  the  face  area  associated  with  the  interface  between  Elements  ( B )  and  (A);  and 
R  is  element  residual,  V  is  estimated  normal  sub-face  flux,  U  is  sub-element-based  quantity 
employed  for  flux  correction,  as  defined  previously.  It  is  noted  that  the  sub-face  area  is  equal  to 
the  face  area  divided  by  the  number  of  face  nodes,  e.g.,  it  is  4  for  quadrilateral  and  3  for  triangu¬ 
lar  faces. 

Scenario  2.  Around  a  boundary  node:  In  Figure  5,  Node  8  is  a  boundary  node  at  which  four 
hexahedral  elements  are  connected.  It  is  associated  with  eight  sub-face  fluxes,  where  four  of 
them  are  associated  with  boundary  faces.  By  setting  Usw  =USX  =USY  =USZ  =  0 ,  the  four  ele¬ 
ment  residual  equations  associated  with  Node  8  are 

p  A  _  rmc  .  rmc  .  rmc 

K8  “  J  8,  BA  +  J8,AW  +  J8,AC 

(2  ^ 

=Av>.‘Mu^-uA}^L+[v^+i(>-uA}^f +[FMc+Kc-f/M)]-J7L 

p B  _  rmc  .  rmc  ,  rmc 

K8  ~  J  8, BA  +  J8,BX  +  J8,BD 

r  /  \  —}  A  r  /  |  A  r  ,  A  (37) 

=  K.e, + K,  ~u>., )]  •  ■ -f- + (°  • -  u%m  )]  ~f-+ [X*, + (u,* -uA}-f- 
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(a) 


(b) 


Element 

(B) 


R“ 


Element 

(D) 


R° 


Figure  5.  The  four  hexahedral  elements  associated  with  a  boundary  node  (i.e.,  Node  8)  and  the  eight 
conservative  sub-face  flows  around  Node  8:  (a)  elements  are  connected,  (b)  elements  are 
separated. 


r>  L)  _  rmc  ,  rmc  r 

K8  —  J  8, DC  J 8, DZ  ~  J 8 


me 
8  ,BD 


^8,DC  (C$,C  ^8,Z>) 


DC 


+ 


^8  ,Z)Z+(b  CgjjJ 


A 


DZ 


v  +(u  -U  'll  •  Abd 

r8,BD^yJi,D  ^  i,B  ) ]  ^ 


(39) 


Because  the  specified  boundary  fluxes  are  already  accounted  for  in  computing  element  residuals, 
Equation  36  through  Equation  39  will  be  modified  if  they  involve  flux-type  boundary  faces.  For 
example,  if  Boundary  Faces  7-8-14-13  and  8-9-15-14  are  assigned  flux-type  boundary  condi¬ 
tions,  Equation  36  and  Equation  37  become 


nA,bf  r>  A  rmc  rmc  ,  r 

n8  ~  K8  J8MW  ~  J8.BA  +  J% 


me 

BA  1  J  8, AC 


^8,&4  "*'(^8,X 

j}B,bf  _  tjB  rmc  _  rmc  .  rn 

K8  ~  K8  ~  J 8, BX  ~~  J 8, BA  +  / 8, 

^8 ,BA  +  (C»,A  _  CS'B  ) 


n  A 


BA 


J  4 


+ 


K ,AC  +  (^8,C  A  ) 


A 


AC 


(40) 


i  A 


BA 


J  4 


+ 


V  +(U  -TJ  )]  ■  ^bd 

r  8,BD  t \^8,D  ^8,B  j 


(41) 


In  this  case,  we  solve  Equation  38  through  Equation  41  for  UgA  ,  US  B  ,  Ug  c ,  Us  D . 


When  all  the  boundary  faces  associated  with  Node  8  are  specified  with  flux-type  boundary  con¬ 
ditions,  Equation  38  and  Equation  39  become 
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TjC,bf  _  nC  rmc  _  rmc  r, 

A8  ~  K8  ~  J$,CY  ~  JS,DC  ~  J 8 


me 
8, AC 


v  +(u  -JJ  ]] .  ^ DC 

K8,Z)C  ^8  ,DJ_\  ^ 


^MC  +  (Cg'C  CS  A  ) 


A 


AC 


r\D,bf  _  t\D  rmc  _  rmc  r, 

K8  ~  K8  _  J 8, DZ  ~  J8,DC  ~  J 8 


me 

8, DC  J8,BD 


y  _| Jtj  —U  ^ DC 

v  8, DC  '  8,C  ^  8,D  ) ]  ^ 


Vg,BD+{Ug'D  Us  B^j 


A 


BD 


(42) 


(43) 


Now  the  four  element  residual  equations  are  Equation  40  through  Equation  43.  But  one  of  them 
has  to  be  replaced  by  setting  one  U  value  to  zero  due  to  linear  dependency  among  the  four  equa¬ 
tions,  as  discussed  before.  If  here  we  set  U&  A  =  0 ,  and  use  it  to  replace  Equation  40,  we  are  to 

solve  Equation  41  through  Equation  43  for  USB,  Usc,  UgD.  The  computed  If  s  are  then  used  to 

correct  nonnal  sub-face  flux,  followed  by  taking  arithmetic  average  of  the  conservative  sub-face 
fluxes  to  obtain  conservative  normal  face  flux. 


The  Global  Approach.  In  the  global  approach,  the  conservative  normal  fluxes  are  first  related 
to  element  residuals  at  the  element  level  to  satisfy  the  local  conservation.  A  global  residual  equa¬ 
tion  (in  matrix  form)  is  then  formed  by  assembling  all  element  residual  equations  and  solved  for 
all  the  element-based  If  s.  These  If  s  are  used  to  correct  the  normal  fluxes  through  element 
edges/faces  for  mass  conservation. 

Scenario  1:  For  elements  whose  edges  are  all  internal.  As  shown  in  Figure  6  below,  Element 
( A )  is  an  internal  element:  all  of  its  element  edges  are  internal.  Therefore,  the  following  element 
residual  equation  exists  for  Element  (A): 


RA=Rf+R2A+R?=fZc 


/'me  .  rmc  j ymc  j 

CA  ~l”  J  AD  ~  V  AB  '  . 


ymc  #  7  ,  ymc  #  i 

v  CA  lCA  '  v  AD  lAD 


(44) 


where 

Re  =  total  element  residual  of  element  (E) 

Rf  =  element  residual  associated  with  Node  i  at  element  (E) 

fyfy  =  normal  edge  flow  from  Element  (Ej)  to  Element  (ZA),  which  conserves  mass  locally 
F~2  =  normal  edge  flux  from  Element  (£/)  to  Element  (ZA),  which  conserves  mass  locally 
lFEi  =  edge  length  associated  with  the  interface  between  Elements  (ZA)  and  (ZA) 


Equation  44  can  be  written  further  as 

«A  =*.,  +  *2,  +  r>a  =[rM  +(UB -ut)\-iM-[va  +(uA -UC)}1CA 
+  [VAD+(UD-U.t)\lAD 


(45) 


where  VE^  =  given  estimated  normal  edge  flux  from  Element  (Ej)  to  Element  (ZA). 
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Figure  6.  An  internal  triangular  element  (i.e.,  Element  (A))  and  its  three  neighbor  elements:  (a)  elements 
are  connected,  (b)  elements  are  separated. 


Equation  45  is  the  residual  equation  for  Element  (A).  Likewise,  we  can  compose  the  residual 
equation  for  each  internal  element.  The  conservative  fluxes  through  the  three  internal  edges  of 
element  (A)  are  computed  with 


Vab+(Ub-Ua) 

(through  Edge  1-3) 

(46a) 

vca+(ua-uc) 

(through  Edge  1-2) 

(46b) 

to 

+ 

c: 

to 

1 

c: 

(through  Edge  2-3) 

(46c) 

Scenario  2:  For  elements  that  contain  boundary  edges.  As  shown  in  Figure  7  below,  Element 
(A)  has  two  boundary  edges  (i.e.,  Edges  1-2  and  1-3)  and  one  internal  edge  (Edge  2-3).  Without 
the  flux-type  boundary  conditions  applied,  the  residual  equation  for  Element  (A)  is: 

D  —  D  i  D  I  D  _  fmc  I  fmc  I  fmc 

‘'a  —  “l/f  +  ^2  A  +  ^3  A  —  Jab  +  Jaz  +  I  AX 

(47) 

+|A  Auz -uA)\iAZ  +|>„+(1/V -UA)\1„ 

Situation  1.  All  boundary  edges  are  assigned  flux-type  boundary  conditions:  When  flux-type 
boundary  conditions  have  been  taken  into  account  in  computing  the  element  residuals,  fff  and 
f'fx  will  no  longer  appear  in  the  element  residual  equations.  If  qAz  and  q,\x  represent  the  outward 
boundary  flux  through  Edges  1-2  and  1-3,  respectively,  Equation  47  becomes: 

RlI  =  Ra  -  -  A  AX  =  fZ  =  [y.,3  +  (U.  -  UA )]  •  (48) 
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Figure  7.  A  boundary  triangular  element  (i.e.,  Element  (A))  that  has  two  boundary 
edges  and  one  internal  edge. 


The  outward  normal  flux  through  Edges  1-2,  1-3  and  2-3  are  computed  by 


Kz  =  qAZ 

(through  Edge  1-2) 

(49a) 

rS  =  qAx 

(through  Edge  1-3) 

(49b) 

Vmc  =  V 

v  AB  v  AB 

+  (UB-UA)  (through  Edge  2-3) 

(49c) 

Situation  2.  Some,  but  not  all,  boundary  edges  are  assigned  flux-type  boundary  conditions:  In  the 
two-dimensional  case  here,  we  consider  when  only  one  of  the  two  boundary  edges  is  assigned  a 
flux-type  boundary  condition.  Suppose  the  outward  boundary  flux  through  Edge  1-2  is  qAz, 
Equation  46  becomes: 

Rbi  =  K-qv  =fZ+/Z 

=  B-u  ,)}iA,+\y„+{Q-u  A)\1AC 

The  outward  normal  fluxes  through  the  three  element  edges  are  computed  with 

V'f  =  qA7  (through  Edge  1-2)  (51a) 
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rs=^+(o  -£/,) 

(through  Edge  1-3) 

(51b) 

Va~  =  Vab+(U,-Ua) 

(through  Edge  2-3) 

(51c) 

Situation  3.  None  of  the  boundary  edges  is  assigned  a  flux-type  boundary  condition: 
the  outward  nonnal  fluxes  can  be  computed  by  using 

In  this  case, 

^=Vaz+(0-Ua) 

(through  Edge  1-2) 

(52a) 

vS  =  rM+(o  -uA) 

(through  Edge  1-3) 

(52b) 

vZ  =  vA,+(u,-uA) 

(through  Edge  2-3) 

(52c) 

Extension  to  three-dimensional  space.  The  extension  of  the  global  approach  from  two-  to 
three-dimensional  space  is  also  straightforward.  Again,  we  use  hexahedral  elements  for  demon¬ 
stration. 


Scenario  1:  For  elements  whose  faces  are  all  internal.  As  shown  in  Figure  8  below,  Element 
(A)  is  an  internal  element:  all  its  six  element  faces  are  internal.  The  element  residual  equation 
associated  with  Element  (A)  is 


R4  =  R4  +  R4  +  R4  +  R4  +  R4  +  R4  +  R4  +  R4 

/'me  .  rmc  .  rmc  .  rmc  rmc  .  rmc 

AB  '  J AC  ^  J  AD  ^  JAB  ^  JAF  '  J AG 

=[v‘b+(u,-  U., )]  •  aab +\yAC+(uc-uA)\  aac + [vAD + (u„  ~uA )]  •  aad 
+  [V^+(UB-UA)\-AAE+[vAF+(UF-Ut)\-AAF+[vAB+(Ua-UA)\-AAB 


(53) 


Scenario  2:  For  elements  that  contain  boundary  faces.  As  shown  in  Figure  9  below,  Faces  2- 
3-7-6,  1-2-6-5,  and  1-4-3-2  that  are  associated  with  Element  (A)  are  boundary  faces.  Without  the 
flux-type  boundary  conditions  applied,  the  residual  equation  for  Element  1  is: 

R4  =  R4  +  R4  +  R4  +  R4  +  R4  +  R4  +  R4  +  R4 

/‘me  ,  rmc  ,  rmc  ,  rmc  ,  rmc  .  rmc 

AB  J  AX  J  AY  J AE  +  J AZ  J  AG 

+[v*+(v*-  UA )]  •  aae  +  [vAZ + (o  ■ -(7, )] .  aaf +[rM+(ua-uA)\-  am 

Situation  1.  All  boundary  faces  are  assigned  flux-type  boundary  conditions:  If  q  ix,  q.4Y,  qAz 
represent  the  specified  outward  boundary  fluxes  through  the  three  boundary  faces  and  have  been 
taken  into  account  in  computing  the  element  residuals,  Equation  54  becomes: 
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Figure  8.  An  internal  hexadedral  element  (i.e.,  Element  (A))  that  has  six  internal  faces. 


Figure  9.  A  boundary  hexadedral  element  (i.e.,  Element  (A))  that  has  three  boundary  and  three  internal 
faces. 
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R‘I  =R*-‘!AX-‘lAr-‘l,z=  fZ  +  fZ  +  fZ 

=  \yA,+(u.-uA)\AM+\yAE+(uE-ui)\-AiE+[rM+(ua-uA)\-AM 
The  nonnal  fluxes  for  Boundary  Faces  2-3-7-6,  1-2-6-5,  and  1-4-3-2  are  computed  with 


vz  =  qAX 

(through  Face  2-2-1 -6) 

(56a) 

Vfy  =  qAY 

(through  Face  1-2-6-5) 

(56b) 

VTz=qAZ 

(through  Face  1-4-3-2) 

(56c) 

Situation  2.  Some,  but  not  all,  boundary  faces  are  assigned  flux-type  boundary  conditions:  Sup¬ 
pose  Face  2-3-1 -6  is  the  only  boundary  face  of  Element  (A)  assigned  an  outward  normal  boun¬ 
dary  flux,  qAX,  the  element  residual  equation  for  Element  (A)  is 

T)A,bf  j}A  „  rmc  ,  rmc  .  rmc  .  rmc  .  rmc 

K  =K  -Vax  =  Jab  +Jad  +Jae  +  J  af  +  Jag 

=  [vM  +  (u,  -£/,)]■  Aab  +  [VjD  +(0-Ua)\aab+[vae+(Ub-Ua)\- Aae  (57) 
+  \VAF  +  (0  -  UA )]  •  Aab  +  [VA0  +  (UG -UA )]  •  AJ0 

The  outward  normal  fluxes  through  Boundary  Faces  2-3-7-6,  1-2-6-5,  and  1-4-3-2  are  computed 
by 

V'ff  =  qAX  (through  Face  2-2-1 -6)  (58a) 


yZ=rA,+(o-uA) 

(through  Face  1-2-6-5) 

(58b) 

yZ=r*+(o-uA) 

(through  Face  1-4-3-2) 

(58c) 

Situation  3.  None  of  the  boundary  faces  is  assigned  a  flux-type  boundary  condition: 
the  outward  nonnal  fluxes  can  be  computed  by  using 

In  this  case, 

rZ=r*+(o-uA) 

(through  Face  2-3 -7-6) 

(59a) 

r:;  =  r,Y+(o-uA) 

(through  Face  1-2-6-5) 

(59b) 

r,7=r,z+(o  -uA) 

(through  Face  1-4-3-2) 

(59c) 

Summary  of  conservative  flux  computation  using  the  global  approach.  We’ve 
discussed  how  the  element  residual  equation  for  an  element,  either  internal  or  containing 
boundary  faces,  is  composed  in  both  two-  and  three-dimensional  spaces.  Suppose  there  are 
totally  M  elements  in  the  FE  computational  mesh;  we  can  construct  M  element  residual 
equations,  one  for  each  element.  These  M  element  residual  equations  are  assembled  to  form  a 
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global  matrix  equation  based  on  the  face-element  connectivity.  This  global  matrix  equation  is 
then  solved  for  the  element-based  IPs  that  are  used  to  correct  normal  face  flux  for  mass 
conservation. 

SOFTWARE  DEVELOPMENT:  We  developed  a  software  toolkit  named  Consistent/inconsistent 
Conservative  Flux  Computation  Toolkit  (CCFlux)  to  assist  application  developers  with  comput¬ 
ing  locally  conserved  fluxes.  Systems,  e.g.,  the  shallow  water  system,  solved  in  the  primitive 
form  would  preserve  local  mass  conservation  when  fluxes  are  computed  consistently.  However, 
systems,  such  as  subsurface  flow,  solved  in  a  derived  fonn  without  holding  an  explicit  conserva¬ 
tion  law  would  not  produce  local-mass  conserved  fluxes  because  the  flux  is  computed  inconsis¬ 
tently  with  the  discrete  equation  being  solved.  The  goal  is  to  obtain  locally  conserved  flux  across 
edges  in  2-D  or  faces  in  3-D  systems.  In  the  CCFlux  toolkit,  functionalities  include  (1)  construc¬ 
tion  of  a  unique  edges/faces  map,  (2)  local  renumbering  of  neighbor  elements  around  a  node  on 
the  index  space,  and  (3)  an  efficient  edge/face  manipulation  for  divergence-free  operation.  The 
toolkit  will  be  incorporated  into  different  application  codes,  e.g.,  ADH  (C  code)  [5], 
pWASH123D  (mixed  C  and  Fortran  code)  [6],  or  even  the  old-fashioned  Fortran-77  FEMATER 
code  [7].  Parallel  mesh  manipulation  for  flux  calculation  is  embedded  in  CCFlux. 

Presumably  the  mesh  has  a  unique  global  element  ID  (GEid)  assigned  to  each  element  shown  in 
Figure  10,  though  each  processor  has  its  own  local  set  of  element  IDs  (LEid).  This  requirement 
defines  unique  flux  direction  systematically  without  exception  in  parallel  simulation.  In  CCFlux, 
the  first  and  second  elements  adjacent  to  an  edge  or  a  face  store  lower  GEid  and  higher  GEid, 
respectively,  but  perhaps  lower  LEid  and  higher  LEid,  respectively.  However,  the  flux  direction 
across  a  boundary  edge  is  always  defined  from  inside  to  outside  of  the  domain.  Therefore,  the 
second  element  adjacent  to  a  boundary  edge  is  “-1,”  implying  nonexistent. 


Figure  10.  A  partitioned  2-D  domain  showing  global  element  ID  and  the  edge  list 
collecting  unique  edges  on  each  processor. 
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In  Figure  11,  the  red  color  arrows  around  vertex  32  show  the  flux  direction  associated  with  each 
adjacent  sub-edge.  The  total  flux  of  Edge  31-32  is  the  sum  of  two  computed  fluxes  of  this  sub¬ 
edge  based  on  node  32  and  node  3 1 .  Note  that  vertex  32  and  vertex  3 1  can  be  owned  by  different 
processors.  The  link  list  sketched  at  the  bottom  of  Figure  1 1  collects  all  the  unique  edges  on  each 
processor.  The  pointer  “Next”  points  to  the  next  unique  edge.  However,  “Next2[0]”  and 
“Next2[l]”  point  to  the  next  edge  adjacent  to  the  respective  vertex  comprising  this  edge.  The 
first  entry  of  “Next2”  points  to  the  next  edge  adjacent  to  the  vertex  with  lower  local  vertex  ID  on 
the  edge,  while  the  second  entry  associated  with  the  other  vertex  on  this  edge. 


Figure  1 1 .  Three  link  lists  collecting  all  unique  edges,  edges  adjacent  to  vertex 
32,  and  vertex  36. 


This  design  of  data  structure  avoids  duplicate  edge  nodes,  which  can  waste  memory  and  CPU 
time  spent  in  copying  data  among  duplicate  edge  nodes.  Each  edge  node  has  two  data  fields  to 
store  computed  data  defined  by  application  codes.  Each  one  is  obtained  based  on  each  vertex  on 
the  edge.  The  CCFlux  toolkit  also  renumbers  the  element  ID  adjacent  to  each  vertex.  For  exam¬ 
ple,  in  Figure  1 1  those  6  elements  adjacent  to  vertex  are  numbered  to  from  0  to  5.  The  toolkit  can 
also  provides  a  set  of  Boolean  values  indicating  whether  the  adjacent  flux  direction  of  a  node 
needs  to  be  flipped  if  the  circulation  condition  is  imposed. 

The  toolkit  equips  a  set  of  API  functions  listed  in  Table  1.  The  following  three  functions  are 
designed  for  different  application  codes,  and  must  be  incorporated  at  the  initialization  stage. 
They  are  adhCCFinitMesh  for  ADH,  pWASHCCFinitMesh  for  pWASH123D,  and 
fwCCFinitMesh  for  FEMWATER.  Three  destroy  functions,  adhCCFdestroyMesh, 
pWASFICCFdestroyMesh,  and  fwCCFdestroyMesh  are  also  required  at  the  finalized  stage. 
The  data  handler,  CCFmesh,  separates  the  flux  calculation  from  others  in  the  application  to 
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make  the  algorithm  implementation  portable.  The  ultimate  goal  of  the  toolkit  development — 
reusable  and  modular  code — can  thus  be  reached. 


Table  1.  API  functions 

API  function 

Return  value* 

face_node=CCF_GET_1ST_VTX_NEIGHBOR_FACE(mesh,i) 

Return  the  1st  face  (i.e.,  face  node  shown  in 
Figure  1 1 )  adjacent  to  vertex  i. 

face_ptr=CCF_GET_FACE_POINTER(face_node) 

Return  the  face  pointer  encapsulating  the 
following  information:  face  number  (i.e.  0,  1, 
or  2  for  a  triangle)  and  the  element  pointer 
that  the  face_node  resides. 

ierr=CCFface_vertices(mesh,face_ptr,vtx_list) 

Return  the  vertex  list  comprising  the  face  that 
the  face_ptr  points  to. 

Sign_flipped=CCF_GET_FACE_SGN_DATA2VTX(mesh,i) 

Indicating  whether  each  flux  direction  across 
adjacent  faces  next  to  node  i  needs  to  be 
flipped  if  circulation  condition  is  imposed. 

data_ptr=CCFget_face_node_data(mesh,face_node) 

Return  the  address  pointing  to  user-defined 
data  associated  with  the  face_node. 

CCFface_neighbor_elements(mesh,face_ptr,elmlD,elmlD+1) 

Return  the  local  element  ID  adjacent  to  the 
face  that  the  face_ptr  points  to.  The  flux 
direction  is  defined  from  elmIDjO]  to  elmlDjlj. 

elm_face=CCFget_face_tuple(mesh,ptr,j) 

Return  the  2  adjacent  compact  element  ID  (0 
to  5  show  in  Figure  1 1 )  of  the  jth  face  next  to 
the  vertex  encapsulated  by  the  pointer  ptr. 

face_node=CCFget_next_vtx_neighbor_face(mesh,face_ptr,face_node,elmlD[0],i) 

Return  the  next  face_node  of  the  current 
face_node  retrieved  from  the  facejist  based 
on  the  info,  of  vertex  i.  Either  the  one  linked 
by  next2[0]  or  next2[1]  is  returned. 

CCFsend(mesh) 

Synchronize  face  data  among  processors 

*  “Face”  is  used  to  represent  2-D  edge  and  3-D  face  in  this  column 

To  demonstrate  the  algorithmic  implementation,  the  local  approach  of  Larson-Niklasson  method 
is  used  to  calculate  the  U  values  described  from  Equation  2  to  Equation  10.  Each  U  value  is 
associated  with  an  element  adjacent  to  the  vertex  shown  in  Figure  2.  For  the  2-D  shallow  water 
flow  in  ADH,  the  estimated  flux  is  the  length  integration  of  hV  (a  face  is  actually  an  edge  in 
2-D),  shown  in  the  first  loop  of  the  pseudo  code  in  Figure  12.  Followed  is  the  loop  visiting  each 
vertex  to  construct  a  local  linear  system,  in  which  the  total  number  of  rows  equals  the  number  of 
elements  including  the  vertex.  The  EFT  direct  solver  is  then  used  to  solve  this  linear  system  with 
unknown  U.  Once  U  is  obtained,  the  mass  conserved  flux,  i.e.,  (hV)mc,  can  be  calculated  based  on 
the  formulation  from  Equation  11a  to  Equation  Ilf.  The  final  flux  value  across  each  face  is 
actually  the  average  of  the  local  conserved  (hV)mc  associated  with  the  vertex  sitting  on  the  edge. 
Before  the  computation  of  the  final  flux,  the  function  CCFsend  must  be  executed  to 
synchronize  the  edge  data  (hV)mc  for  summing  on  the  node  star.  To  verify  the  result,  we  can 
simply  compute  the  difference  between  the  total  flux  across  all  the  faces  and  the  sum  of 
boundary  fluxes  described  in  the  governing  equation.  This  verification  function  checks  every 
element  throughout  the  entire  domain.  The  maximum  value  of  the  difference  is  then  found  and 
given  to  the  user.  Ideally  it  is  zero  but  can  practically  be  near  the  user’s  prescribed  error 
tolerance  for  the  application  simulation. 
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Foreach  face  do 
data_ptr [0] =n*jN [0] (hV)dl 
data_ptr [ 1 ] =n* Jn [ 1 ] (hV)dl 
Endfor 

Foreach  owned  vertex  do 

pick  elmID[0]  or  elmID[l]  based  on  sign_f lipped 
compute  element  residual,  i.e.,  RXA, 
construct  A  and  b  =>  AU  =  b 

set  U1;A=0,  LU  direct  solver  is  used  to  obtain  sol. 
compute  (hV)mc  for  each  face 
Endfor 

CCFsend (mesh)  /*  sync  face  data  */ 

Figure  12.  Pseudo  code  for  the  Larson-Niklasson  implementation. 

The  implementation  on  a  3-D  mesh  should  be  similar  except  that  some  geometric  issues  need  to 
be  tackled.  The  CCFlux  toolkit  needs  to  be  flawless  for  both  2-D  and  3-D  meshes.  Currently,  we 
are  in  the  middle  of  debugging  3-D  cases.  A  new  version  of  CCFlux  will  be  released  after  the 
3-D  bug-free  implementation  is  completed. 

EXPERIMENTAL  RESULTS:  We  demonstrate  the  locally  mass-conservative  flux  calculation 
on  an  example  used  in  the  article  by  Tate  et  al.  [8].  The  domain  of  interest  is  San  Diego  Bay  in 
CA,  which  covers  an  area  of  approximately  560  km“  and  has  a  mean  depth  of  6.5  m.  This  area 
was  discretized  with  10,999  triangular  elements  and  6,311  nodes,  as  shown  in  Figure  13.  The 
simulation  took  the  tide  shown  in  Figure  14  as  the  driving  force  and  used  a  time-step  size  of 
360  sec.  The  total  simulation  time  was  24  hours,  which  is  about  a  diurnal  cycle.  The  tidal  boun¬ 
dary  condition  was  applied  to  the  nodes  included  in  the  section  of  open  boundary  (Figure  13). 
The  rest  of  the  domain  boundary  was  set  to  closed  boundary,  i.e.,  no  normal  flow. 


Open  Boundary  (Tidal  BC) 

Pacific  Ocean 


6,311  nodes 
10,999  elements 
3  materials:  channel  (orange) 

outside  channel  (brown) 
ocean  (blue) 


Figure  13.  Domain  discretization  of  San  Diego  Harbor. 
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We  selected  six  nodes  to  examine  the  computed  conservative  fluxes.  Among  these  six  nodes, 
Nodes  21  and  22  are  on  the  open  boundary  (Figure  15a),  Nodes  4684  and  4685  are  on  the  closed 
boundary  (Figure  15b),  and  Nodes  405  and  2997  are  interior  nodes  (Figure  15c).  Figure  16 
shows  the  flow  direction,  which  is  defined  from  the  lower  global  element  ID  to  the  higher  one. 
Figure  17  shows  the  water  depth  color  contour  and  scaled  velocity  vector  at  simulation  time  = 
3,600  s.  Table  2  shows  the  estimated  nonnal  sub-edge  flow  (column  2),  the  conservative  nonnal 
sub-edge  flow  (column  3)  and  the  total  nodal  flow  (column  4)  at  that  time  step.  The  flow  rate 
shown  in  columns  2  and  3  are  positive  when  the  nonnal  flow  direction  is  the  same  as  specified  in 
Figure  16.  The  positive  and  negative  signs  associated  with  the  total  nodal  flow  in  the  “Nodal 
Flow”  column  represent  the  out-going  and  in-coming,  respectively,  flow  at  the  nodes.  As  shown 
in  Table  2,  at  simulation  time  =  3,600  the  total  nodal  flow  at  Nodes  21  and  22  are  negative,  indi¬ 
cating  incoming  flow  at  those  two  open  boundary  nodes  due  to  the  rising  tide  (Figure  14).  On  the 
other  hand,  the  total  nodal  flow  at  both  the  two  interior  nodes  (i.e.,  Nodes  405  and  2997)  and  the 
two  closed  boundary  nodes  (i.e.,  Node  4684  and  4685)  was  zero,  as  expected  in  our  derivation 
given  in  this  technical  note.  Therefore,  we  have  verified  the  implementation  of  the  local 
approach  for  computing  conservative  nonnal  fluxes  on  element  edges  in  solving  2-D  shallow 
water  flow  in  ADH.  The  computational  results  show  the  maximum  difference  computed  by  the 
verification  function  mentioned  in  the  section  of  “SOFTWARE  DEVELOPMENT.” 
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(a) 


(b) 


Figure  1 5.  Locations  of  the  six  nodes  selected  to  examine  the  computation  of  conservative  flux: 

(a)  Nodes  21  and  22  are  on  the  open  boundary,  (b)  Nodes  4684  and  4685  are  on  the  closed 
boundary,  and  (c)  Nodes  405  and  2997  are  interior  nodes. 
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Figure  16.  Zoom-in  of  the  mesh  discretization  abound  the  six  selected  nodes:  (a)  Nodes  21  and  22, 

(b)  Nodes  4684  and  4685,  (c)  Node  405,  and  (d)  Node  2997;  numbers  in  red  are  the  global 
node  IDs;  black  arrows  are  used  to  define  positive  flow  direction  across  each  sub-edge 
around  the  selected  nodes. 
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Water  Depth  [m]  :  3600.0 


Figure  17.  Computed  water  depth  and  velocity  vector  at  simulation  time  =  3,600  s. 
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Table  2.  Results  of  sub-edge  and  nodal  flow  associated  with  the  six  selected 
nodes  at  simulation  time  =  3,600  seconds 

Estimated  Flow*  (ft3/s)  Conservative  Flow**  (ft3/s) 

Around  Node  21: 

2.7838E+04 
1.3247E+04 
3.3049E+03 
-1.5188E+04 
-2.1510E+04 

Around  Node  22: 

3.4450E+04 
2.0206E+04 
2.481 7E+03 
-1.5343E+04 
-2.8987E+04 

Around  Node  405: 

4.8676E+02 
6.1 192E+02 
4.3372E+02 
-1.4027E+02 
-7.7235E+02 
-5.2722E+02 

Around  Node  2997: 

-3.3387E-02 
-1.2641E-01 
3.6237E-02 
1.3369E-01 
4.8901  E-02 

Around  Node  4684: 

-4.3777E-06 
-2.6924E-04 
4.41 10E-03 
-7.0438E-03 

Around  Node  4685: 

7.6636E-03 
3.7038E-03 

*  Estimated  flow  =  (computed  water  depth)  x  (computed  velocity)  x  (sub-edge  length) 

**  Conservative  flow  =  (computed  conservative  flux)  x  (sub-edge  length)  using  the  local  approach 

***  Nodal  flow  =  computational  result  by  substituting  the  convergent  solution  into  the  matrix  equation  of  mass  conservation 


SUMMARY:  In  this  technical  note,  two  algorithms  to  compute  locally  mass-conservative 
edge/face  fluxes  are  presented.  The  local  approach  wins  over  the  global  approach  and  these  were 
derived  based  on  the  Larson-Niklasson  method.  Obviously  the  global  approach  requires  more 
memory  consumption  for  the  linear  system,  contrasted  to  the  local  approach,  whose  linear  sys¬ 
tem  size  is  defined  by  the  total  number  of  adjacent  elements  to  a  node.  The  global  approach  may 
require  an  efficient  linear  solver,  while  the  local  approach  can  simply  use  a  direct  solver. 
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Our  software  development  aims  to  develop  the  CCFlux  toolkit  to  support  the  solution  of  the  local 
approach.  A  set  of  light-weight  application  interface  functions  facilitates  easy  incorporation  to 
different  applications.  The  aforementioned  local  approach  was  implemented  in  the  2-D  shallow 
water  flow  module  in  the  ADH  model.  An  experimental  case  solving  the  San-Diego  Bay  area 
was  presented  for  demonstration.  The  result  verified  the  correct  implementation  of  the  locally- 
conservative  flux  computation  and  the  CCFlux  toolkit.  More  cases,  including  3-D  groundwater 
flow,  3-D  CFD  problems,  and  3-D  shallow  water  flow,  will  be  built  for  demonstration  in  the  near 
future. 

ADDITIONAL  INFORMATION:  This  work  was  funded  by  the  SWWRP.  At  the  time  of  publi¬ 
cation,  the  SWWRP  Program  Manager  was  Dr.  Steven  L.  Ashby:  Steven. L.Ashby@usace. army, 
mil.  This  technical  note  should  be  cited  as  follows: 

Cheng,  J-R.  C.,  H-P.  Cheng,  M.  W.  Farthing,  and  C.  E.  Kees.  2010.  Computing 
locally-mass-conservative  fluxes  from  multi- dimensional  finite  element  flow 
simulations.  ERDC  TN-SWWRP-1 0-4.  Vicksburg,  MS:  U.S.  Army  Engineer 
Research  and  Development  Center,  https://swwrp.usace.army.mil/. 
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